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1.2  Reserach  Objectives  Which  Have  Been  Accomplished 

A  copy  of  the  above-mentioned  report  is  included  with  this  final 

A 

report.  A  paper  entitled  "Curved  Space  Scattering"  has  already  been  sent 
to  AFOSR.  The  abstracts  of  the  paper  and  report  contain  statements  of 
the  research  accomplished.  In  addition  we  have  looked  into  the  following 
problems: 

<— ■  a.  We  have  employed  the  Kanal-Moses  variational  principle  to  treat 
the  synthetic  data  discussed  in  the  report  listed  in  1.1.  It  was  found 
that  if  the  initial  trial  function  is  within  10%  of  the  actual  result, 
then  the  K-M  variational  principle  gives  the  result  to  better  than  1%. 

b.  We  have  generalized  the  results  of  Kay  for  n-poles  in  such  a 
way  that  practical  applications  are  possible,  e.g.  to  the  ionosphere.  To 
test  the  method  we  are  in  the  process  of  studying  the  4-  and  10-pole 
cases  before  treating  the  100-pole  case. 

/  _  -  ,T  ' 
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1.3  Interaction  with  Other  Investigators 

a.  We  discussed  scattering  at  a  Physics  Colloquium,  University 
of  Pennsylvania,  February,  1980. 

b.  We  have  been  discussing  problems  of  mutual  interest  on  a 
continuing  basis  with  Prof.  H.  Moses  of  the  University  of 
Lowell  and  Dr.  A.  Skalafuris  of  the  Naval  Research  Laboratory. 
The  problem  of  obtaining  information  beyond  a  potential  peak 
was  discussed  with  Prof.  Percy  Deift  of  the  Courant  Institute 
and  with  Dr.  Robert  Greene  of  Science  Applications  Inc. 

We  have  had  extensive  discussions  concerning  inverse  scat¬ 
tering  with  Prof.  C.  V.  Vishveshwara  of  the  Raman  Institute 
in  Bangelore,  India.  A  collaboration  between  India  and  the 
United  States  is  now  under  way. 

Applications  to  the  ionosphere  have  been  discussed  with  Drs. 

A.  Jordan  and  S.  Ahn  of  the  Naval  Research  Laboratory.  Appli¬ 
cations  to  oceanography  have  been  discussed  with  Dr.  E.  Toton 
of  the  Naval  Oceanographic  Laboratory  and  Dr.  R.  Adams  of  the 
Applied  Physics  Laboratory. 


A  COMPARISON  OF  THE  APPROXIMATE  AND  EXACT  FULL  WAVE  THEORY 
FOR  THE  SOUNDING  OF  A  STRATIFIED  IONOSPHERE 


Jeffrey  Cohen  and  Michael  Kearney 
Department  of  Physics#  University  of  Pennsylvania 
Philadelphia#  Pennsylvania  19104 


1.  INTRODUCTION  AND  SUMMARY 

The  purpose  of  the  present  communication  is  to  compare 
methods  of  obtaining  the  electron  density  profile  N(z),  which  is 
assumed  to  be  a  function  of  altitude  z,  using  the  usual  approxi¬ 
mate  W.K.B.  method  and  using  the  full-wave  method.  The  data 
needed  to  compute  N(z)  using  the  two  methods  is  identical.  This 
data  is  the  time  as  a  function  of  frequency  which  it  takes  a 
train  of  horizontally  polarized  electromagnetic  waves,  trans¬ 
mitted  vertically#  to  be  reflected  from  the  ionosphere  back  to 
the  transmitter. 

Since  the  data  used  is  identical  for  the  approximate  and 
full-wave  theories,  there  is,  in  principle  at  least#  no  need  to 
modify  equipment.  The  difference  in  treatments  is  essentially 
computational  and  the  same  experiment  can  be  used  to  obtain 
N ( z)  using  the  approximate  and  exact  theories. 

The  approximate  treatment  has  been  in  use  for  over  forty 
years  and  many  of  the  purely  computational  difficulties  have 
long  since  been  overcome.  Those  who  use  ionosondes  and  calculate 
profiles  have  become  so  accustomed  to  the  use  of  the  approximate 
method  that  they  often  forget  the  method  really  is  approximate 
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The  exact  treatment  suffers  from  the  defect  that  it  is  new  and 
therefore  the  purely  computational  methods  are  only  now  being  set 
up.  Indeed  one  of  the  objects  of  the  present  communication  is  to 
stress  the  need  for  understanding  the  computational  difficulties 
and  thus  make  it  possible  to  apply  the  full  wave  theory  to  the 
data.  But  the  principal  point  of  this  introduction,  and  of  this 
communication  is  to  provide  a  larger  view  of  the  problem. 
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2.  THE  EXACT  DIFFERENTIAL  EQUATION  FOR  THE  ELECTROMAGNETIC 
FIELD. 

Let  E(Z/t)  be  a  component  of  the  electric  field  perpen¬ 
dicular  to  the  z-axis  or  vertical  axis  along  which  the  electron 
density  is  stratified.  Then  with 
Ek(z,t}  =  elwt  Ek(z)  , 

(w  =  ck)  (1) 

the  equation  for  E^(z,t)  is  (ref.  1,  page  129) 

—  E.  (z)  +  k2  n2  E,  (z)  =  0  (2) 

dz2  *  k 

In  (1)  a)  is  the  angular  frequency,  k  =  —  is  the  wave  number  in 

c 

free  space  and  n  is  the  index  of  refraction  given  by 
n2  =*  1  -  N(z)  4tte2 

mc2k2  (3) 

where  m,  e,  c  are  the  mass  of  the  electron,  the  charge  of  the 
electron  and  velocity  of  light  in  Gaussian  cgs  units. 

In  our  discussion  we  have  ignored  electron  collisions  and 
the  effect  of  the  earth's  magnetic  field.  This  physical  (as 
opposed  to  mathematical)  approximation  is  often  assumed  in 
ionospheric  sounding  methods.  At  the  magnetic  equator,  the  earth's 
magnetic  field  is'"  parallel  to  the  earth's  surface  and  hence  the 
polarization  of  the  ionosonde  can  be  chosen  so  that  the  magnetic  field 
does  not  affect  the  scattered  wave.  In  other  situations,  the  effect 
of  the  magnetic  field  may  also  be  eliminated  (see  Ref.  3) .  The 
effect  of  collisions  is  negligible  at  higher  altitudes  (z  >  80  km) . 

One  of  the  ultimate  objectives  of  the  present  train  of  research  is  to 
determine  whether  these  physical  approximations  lead  to  errors  greater 
than  those  made  by  the  mathematical  approximations.  In  any  case, 
we  shall  assume  as 
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is  customary  that  the  use  of  Eq.  (3)  for  the  index  of  refraction 
is  a  good  approximation  for  use  in  ionospheric  sounding  in  certain 
geographical  regions. 

Equation  (2)  can  also  be  written 

— —  (z)  +  (k2  -  V(z) )  Ek(z)  =  0  (4) 

dz2 

where  V(z)  is  given  by 

V(z)  =  K  N(z),  with  K  =  . 

me2  (4a) 

Equation  (4)  is  the  one-dimensional  Schroedinger  equation 
which  has  been  exhaustively  studied.  There  has  been  a  renewal  of  inter' 
arising  in  recent  years  because  of  its  connection  with  soliton 
theory  (4) .  The  potential  V(z)  is  essentially  the  number  density 
N(z)  . 

In  the  direct  problem  of  reflecting  electromagnetic  waves 
from  the  ionosphere  we  assume  V(z)  (or  N(z))  is  very  small  for 
z  <  Zq  and  z  >  z^. 

We  look  for  solutions  of  (4)  which  behave  like 

Ek(z)  =  elkz  +  b (k)  e"lkz  for  z  <  zQ 

=  t(k)  e^kz  for  z  >  z^  (5) 

The  quantities  b(k)  and  t(k)  are  called  the  reflection 
and  transmission  coefficient  respectively. 

For  a  wave  with  wave  number  k, 
we  have 

EU,t)  =  eik(z-ct)  +  b(k)  e-ik(z+ct)i  z  ^ 

.  .  ik(z-ct)  _ 

=  t (k)  e  ,  z  >  z^ .  (6) 
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The  boundary  conditions  (b)  are  interpreted  to  mean  that 

2TT 

a  plane  wa'.e  with  wave  number  k  =  -j-  moves  initially  toward 
the  scattering  potential  V(z)  and  then  is  partially  reflected 
toward  the  left.  A  transmitted  wave  on  the  other  side  of  the 
potential  which  moves  to  the  right  is  also  present. 

Actually#  one  never  sends  an  infinite  plane  wave  toward 
the  potential  (or  ionosphere) .  Instead  one  sends  a  pulse  con¬ 
taining  several  wave  lengths.  This  pulse  can  be  represented  by 

+°° 

E(X/t)  =  /  A(k)  e"iaJt  Ek(z)  dk  (7) 


and  is  thus  a  superposition  of  the  infinite  plane  waves.  The 
amplitude  factor  A(k)  has  its  peak  value  near  or  at  the  value 
of  k  =  kQ ,  where  kQ  is  the  wave  number  of  the  plane  wave  which 
appears  in  the  pulse.  One  sends  in  a  pulse 

incident' <*•«  =  /"»»)  eik(x-ct)  dk  (8) 

— oo 

and  gets  back  a  reflected  pulse 

Erefl«cted(x't)  '  /”*<*>  »(k)  e-ik(x+ct)  dk.  (9) 

_oo 

since  Eincident(x't)  and  Reflected (x,t)  are  of  finite  extent 
we  can  ask  for  the  time  for  the  reflected  wave  to  return  to  the 

transmitter.  This  time  will  depend  on  kg.  It  is  given  by 

T(kg)  =  |  evaluated  at  k  =  kQ  (10) 

where  <j»(k)  is  the  phase  of  b(k);  i.e. 


b (k)  =  b  (k)  I  e 


i<Mk) 
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The  virtual  height  h^  is  given  by 

hv  =  |cT(kQ).  (11a) 

It  should  be  remembered  that  the  theory  thus  far  is  exact. 

(In  using  (10)  for  the  time  delay  it  is  convenient  to  think  of 
the  transmitter  as  being  located  at  z  =  0.) 
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3.  THE  APPROXIMATE  INVERSE  PROBLEM 

The  direct  problem  is,  for  our  purposes ,  the  following: 
Given  V(z)  (or  equivalently  N(z))/  find  T(k)  . 

The  inverse  problem  is:  Given  T(k)  for  all  k,  find  V(z) 
(or  N(z) )  . 

The  approximate  solution  is  obtained  from  the  W.K.B.  ap¬ 
proximation.  Assume  ( z )  has  the  form 

Ek(z)  =  A  elF(z)  .  (12) 

Substitute  (12)  into  equation  (2)  to  obtain 


,dF.  2 


<■&  m  k2  n2  +  1 


.  d2F 


which  becomes  a  Riccatfc  equation  if  y  =  dF/dz.  An  iteration  pro¬ 
cedure  can  be  started  by  assuming  on  the  right  of  eq.  (13)  that 
2  2  d2F 

k  n  is  large  compared  to  -  ,  i.e.  n/X  is  large  or  n  varies 

dz2 

slowly  compared  to  the  wave  length  (see  Ref.  3  for  a  more  careful 
analysis  of  domain  of  validity) .  The  first  iteration  gives 
F'  =  +  kn  while  the  second  yields 

•It  It  t 

ip  *5  i  p  in 

F'  =  +  kn  (1  +  — - )  ~  +  kn(l  +  — - )  =  +  kn  + 

k2n2  “  2k2n2  "  2n 

which  gives  F  by  a  quadrature.  Since  Eq.  (13)  is  non-linear/ 

the  two  particular  solutions  of  (13)  cannot  be  superimposed  to 

find  a  general  solution.  However,  the  solution  F  can  be  substituted 

into  Eq.  (12)  to  obtain  particular  solutions 

Efc  ( z)  =  A[n(z)J  **.  exp  ["+  i  k  /  n(z')  dz'J  .  (14) 

z0 

to  Eq.  (2).  If  one  particular  solution  to  Eq .  (2)  is  known,  it 
can  be  used  to  give  a  first  order  linear  equation  which  can  be 
solved  by  quadratures  thereby  giving  the  general  solution. 
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Alternatively,  two  distinct  particular  solutions  can  be  added 
to  give  the  general  solution  of  Eq.  (2)  which  is  (within  this 
approximation) 

E.  (z)  =  A(k)  £n(zl]  exp  [+  i  k  /  n(z')  dz'J. 

2° 

+  B  (k)  [n(z)J  1//2  exp  [-  i  k  f  n(z')  dz'J. 

z.  (■ 

Using  the  first  boundary  condition  of  Eq.  (6)  and  noting 
n(z)  =  1  for  z  ^  zQ,  one  finds 
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Ek(z)  =  elkz  +  b(k)  e"lkz  (z  <  zQ) 

=  A(k)  e”ikz0  eikz  +  B (k)  eikZ°  e'ikz 


ikz, 


Hence  A(k)  =  e 


h(k)  =  b(k)  e 


ikz. 


(16) 


Let  z^  be  the  value  of  z  for  which  n(z)  =  0,  i.e.  from  Eq.  (3) 


and  (4a) 


V(zk)  =  k2  . 


(17) 


Our  picture,  the  usual  one,  is  that  V(z)  =  0  for  z  <  zQ  and  • 
increases  monotonically  toward  a  maximum  as  z  increases  above 
Zq.  As  k2  varies  from  0  to  the  maximum  of  V(z^),  a  z^  is  defined 
by  (17). 

From  (15),  Ek(zk)-*  «*>  unless 

2k  -  -  z’ 

r 

2, 


A (k)  exp[i  k  /  n(z’)dz'J  +  B(k)  exp  i  k  /  n(z')  dz'J 

Z0 

(18) 


=  0. 


Thus 


B(k)  =  -  A(k)  expfeik  /  n(z')  dz'J 


ikzn  r  Zk  1 

=  -e  exp[2ik  /  n(z')  dz'J  . 

Or  from  (16)  Z° 

2ikz  p  zk  “T 

b(k)  =  -  e  0  exp]_2ik  /  n(z')  dz'J  . 

20 

The  phase  <})(k)  of  Eq.  (11)  is 


<f>(k)  =7 r  +  2  k 


r  k  -j 

[z  +  /  n(z')  dz'J  . 


(19) 


(20) 


(21) 
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T(k) 


c  *.z0  +  /  nCz'Jdz* 
Z0 

.  2k  ,  ■  8zk 

+  c”  n(zk>  W~  ‘ 


J  +  zr  f  |k  n<z'>  dz' 


(22) 


But  n(zk)  =  0  by  definition  (see  Eqs .  (3)  and  (4)),  also 


3n 


{  -t  =  1_  -L  /v- 2  -v  ( z) 
ak  3k  k  ’'**  vlz; 


=  -  —  /k2-V(z)  +  ---- - 

k2  /k2-V(z) 


=  -  En(z)  +  EhW 


(23) 


Thus  finally 


T(k)  =  §  z  +  / 


dz 1 


n(z') 


or 


4 

(k)  -  f  [«a  ♦  k  / 


zo  ^ 


— 7  . 

-V(z') 


(24) 


This  is  a  well-known  integral  equation  for  V(z) .  The  first 


term  on  the  right 


2z, 


represents  the  time  for  the  signal  to 


reach  Zq  and  return  to  the  transmitter  at  z  =  0.  If  k^  is 
the  lowest  value  of  k  for  which  there  is  a  reflection,  then 


Z0  = 


cT  (kQ) 


(25) 


Now  that  Zq  is  determined,  one  can  find  V(z)  (or  N(z))  using 
Abel's  method  of  inversion  (Ref.  5).  For  each  value  of  k  one 


gets  T(k)  experimentally.  One  then  can  find  V(z)  (zQ  <  z  <  z^) 


using  Abel's  method.  This  is  the  usual  procedure 
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Note  if  k2  >  V  „  (z) /  the  method  fails.  Even  if  k2  is 
max 

near  but  less  than  V  .  we  are  outside  the  domain  of  validity 

max  ■* 

because  V(z)  changes  rapidly  near  its  maximum. 
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4.  THE  EXACT  INVERSE  PROBLEM 

The  exact  or  full-wave  theory  version  of  the  inverse  prob¬ 
lem  is  developed  in  Refs.  1  and  2. 

Let  b(k)  be  the  reflection  coefficient  as  in  Eq.  (6) .  The 
potential  can  be  recovered  from  a  knowledge  of  b(k)  through  the 
use  of  the  Gelfand-Levitan  algorithm.  To  be  specific  let  us 
define  R(z)  by 

*4*00  • 

R(z)  =  (2ir)_1  /  b(k)  e"lkz  dk.  (26) 

—00 

The  reflection  coefficient  b(k)  satisfies  the  following  con¬ 
ditions  (see  Refs.  1,  2  and  6.) 

b(-k)  =  b*(k)#  b(k)  analytic  in  the  upper  half-plane , 

2ikzn 

b (0)  =  -1,  R(z)  =  0  for  z  <  2zq,  b(k)  =  e  g(k).  (27) 

The  phase  of  g (k)  approaches  a  finite  limit  as  |k|-*-®  on  the  real  axis. 

The  time  delay  obtained  from  the  phase  of  b(k),  approaches  the 
2z0 

value  — —  as  |k|  •*■<*>,  i.e.  if  <J>(k)  is  the  phase  of  b(k)  as  in 
c 

Eq.  (11) 

4>(k)  -*■  2kz^  as  |k|  -*■  <*>  .  (28) 

Let  us  define  the  Gelfand-Levitan  kernel  K(z,y)  by 
K(z,y)  =  0  for  either  y  >  z,  or  z  <  z^.  (29) 

For  y  <  z  and  at  the  same  time  z  >  we  require  K(z,y)  to 
satisfy  the  Gelfand-Levitan  equation: 

x 

K(z,y)  =  -  R(z+y)  -  n(z+y-2z  )  /  K(z,u)  R(u+y)  du  (30) 

2*a-y 

where  n(x)  is  the  Heaviside  function  defined  by  n(x)  =  1  for 
x  >  0,  n(x)  =  0  for  x  <  0. 


Having  found  the  Gelfand-Levitan  kernel  K(z,y),  V(z)  is 
given  by  the  simple  expression 

V(z)  =  2  K ( z , z)  .  (31) 

However#  the  electric  field  E^(z)  can  also  be  obtained  using 

E  (z)  «  elkz  +  b (k)  e"lkz  +  j  K(z,u)  [elku+  b(k)  e“lkuJ  du. 

2z  Q-z 

(32) 

Thus  to  find  V(z)  using  the  algorithm,  we  obtain  the  re¬ 
flection  coefficient  b(k)  from  its  phase  <J>(k). 

This  is  accomplished  using  the  analytic  properties  of  b(k) 
and  hence  log  b(k)  and  a  generalized  form  of  the  Hilbert  trans¬ 
form.  To  be  specific,  (d/dk)<j>(k)  is  found  from  the  time  delay 
in  accordance  with  Eq.  (10).  From  Eq.  (28)  z^  can  be  found. 

The  phase  <j>(k)  is  obtained  by  integration  with  the  boundary 
condition  4>(°)  =  it.  Let  v(k)  be  defined  by 

v (k)  =  <j>(k)  -  2kzQ  (33) 

and  w(k)  by 

w(k)  =  log  |b(k)  |  ^  (34) 

then 

(35) 

■ago 

In  Eq.  (35)/  the  symbol  P  means  the  principal  part  of  the 
integral.. 

It  should  be  mentioned  that  there  are  variational  prin¬ 
ciples  available  (Ref.  7)  which  enable  one  to  obtain  V(z)  when 
b(k)  is  known.  These  principles  have  an  upper  bound  built 
into  them. 
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One  also  has  available  a  generalization  of  the  Gelfand- 
Levitan  algorithm.  If  for  a  given  reflection  coefficient 
bg(k)  one  knows  the  electron  density  VQ(z),  one  can  obtain 
V(z)  -  VQ ( z)  in  terms  of  b(k)  -  bQ(k)  (Ref.  8) .  One  can  view 
this  generalization  as  offering  at  least  two  options.  One  may 
regard  V(z)  -  Vg(z)  as  the  error  in  density  due  to  an  error 
b(k)  -  hg(k)  in  the  reflection  coefficient.  Or  one  may  think 
of  Vg (x)  as  being  the  density  associated  with  a  time  delay 
leading  to  b^ (k)  having  been  obtained  from  a  model  or  a  pre¬ 
vious  calculation  (even  using  the  WKB  method) .  Then  V(x)  is 
obtained  as  a  relatively  small  change  due  to  the  change  in  the 
reflection  coefficient.  The  variational  principle  can  also  be 
used  to  obtain  V(x)  -  VQ(x)  from  b(k)  -  bQ (k)  together  with  a 
bound  on  this  difference.  We  refrain  from  details. 
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5.  NUMERICAL  SOLUTION  OF  THE  EXACT  INVERSE  PROBLEM 

A  straightforward  technique  for  solving  the  Gelf and-Levitan 
integral  equation  numerically  is  to  replace  equation  (30)  with  a 
discrete  approximation  on  a  uniform  coordinate  grid  in  the  (z,y) 
plane.  (Refs.  10  and  11) .  Assuming  nonvanishing  z^  and  grid 


size  A  yields  the  approximation 


m+n 

E 

i=0 


K  =  -  R  ,  -  A  E  K  .  R.W.  ,  -m<n<n 
mn  m+n  ,*_rt  m,i-n  l  i,m+n 


(36) 


where  K  =  K(mA,nA)  and  R  =  R(mA) .  The  term  W.  is  a 

mn  m  i  ,m+n 

weighting  factor  which  depends  on  which  quadrature  rule  is  chosen 
to  approximate  the  integral  in  equation  (30) .  The  best  choices 
would  be  those  forms  which  maximize  the  precision  of  the  approxi¬ 
mation  for  a  given  grid  size  and  integrand  and  which  allow  stable 
numerical  computation.  Note  that  the  weighting  factor  may  depend 
on  the  number  of  points  summed  over.  For  example  many  quadrature 
rules  depend  on  whether  the  number  of  grid  points  summed  over  is 
odd  or  even. 

The  discrete  fourier  transform  (Ref.  12)  provides  a 

natural  framework  for  calculating  the  values  of  R  ,  the  fourier 

m 

transform  of  the  reflection  coefficient,  b(k),  at  the  coordinate 
grid  points.  This  follows  from  the  fact  that  fourier  integrals 
can  be  approximated  by  discrete  sums  over  a  uniform  grid  and  that 
computationally  efficient  algorithm,  the  fast  fourier  transform 
(FFT) ,  (Ref.  14)  exists  for  evaluating  these  sums  (Ref.  13). 

Once  the  values  of  R^ji-O,  M  have  been  calculated  (we  assume 
here  that  b(k)  is  given  or  measured  experimentally)  Eq.  (36)  is 


16 


equivalent  to  a  set  of  2m+l  simultaneous  equations  for  each 

value  of  m.  Solution  of  these  equations  will  yield  all  the 

values  of  K  for  -m  £  n  £  m.  If  desired,  scattering  potential 
mn 

may  then  be  derived  from  tie  values  of  K  ,  using  eq.  (31). 

mm 

There  are  a  variety  of  techniques  available  for  solving  such 

sets  of  simultaneous  equations.  These  techniques  may  generally  be 

divided  into  two  classes,  direct  and  iterative.  Direct  methods 

have  the  advantage  of  assured  solution  for  nonsingular  matrices 

but  require  ~N  arithmetic  operations  for  their  solution  (for 

matrices  of  size  N) .  This  is  a  distinct  disadvantage  for  large 

sets  of  equations.  Iterative  techniques  (Gauss-Seidel ,  Jacobi, 

2 

for  example) ,  generally  require  KN  operations  (k  is  the  number 
of  iterations)  but  do  not  assure  convergence  in  less  than  N 
iterations  and  so  may  not  be  superior  to  direct  methods . 
Accelerated  convergence  might  be  accomplished  through  the  use  of 
relaxation  methods,  a  common  approach  in  the  numerical  solution 
of  partial  differential  equations.  However,  the  best  technique 
probably  depends  to  a  great  extent  on  the  problem  at  hand. 

The  authors  have  used  both  approaches  successfully  on  the 
model  calculations  of  the  next  section. 
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6.  NUMERICAL  SOLUTION  FOR  SOME  MODEL  POTENTIALS 

Our  approach  to  testing  the  numerical  methods  described  in 
the  previous  section  is  to  use  synthetic  data  obtained  from 
simple  potential  models  for  which  exact ,  analytic  solutions 
are  available  for  the  potential,  reflection  coefficients  and 
eigenfunctions.  This  has  the  advantage  that  the  potential  V(x) 
and  hence  the  Gelfand-Levitan  kernel  K(x,x)  can  be  evaluated 
easily.  This  gives  a  direct,  simple  technique  for  evaluating 
numerical  methods  for  accuracy,  stability  and  speed. 

We  have  chosen  to  start  with  model  reflection  coefficients 

derived  from  known  potentials  rather  than  from  arbitrary  reflection 

coefficients  corresponding  to  a  potential  with  properties  not 

known  a  priori.  The  solutions  of  the  Gelfand-Levitan  integral 

equation  may  then  be  compared  with  the  known  solutions  in  a 

straightforward  way. 

in  table  I 

Below/,  we  list  the  models  we  have  used  for  preliminary  tests 
of  our  numerical  techniques.  We  give  both  the  potential  form 
and  reflection  coefficient.  For  each  potential  we  have  calculated 
the  virtual  height  (Eq.  (11a)  and  the  fourier  transform  of  the 
reflection  coefficient  for  parameters  which  are  typical  of  those 
observed  in  the  ionosphere:  penetration  frequencies  of  the  order 
~  3mhz,  widths  of  order  -50  km.  These  results  are  presented  in 
Figs.  1  through  14.  One  feature  of  the  fourier  transforms  should 
be  mentioned.  For  the  more  realistic  ionospheric  models  n,  the 
rectangular  and  the  parabolic  electron  distributions,  the  trans¬ 
forms  exhibit  rapid  oscillations  whose  wavelength  is  roughly  Kp“^, 
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where  Kp  is  the  penetration  wave  number.  This  feature  will 
place  a  severe  strain  on  calculations  using  the  discrete 
Gelfand-Levitan  equation,  as  the  grid  size  must  be  made  suffi¬ 
ciently  small  so  that  the  approximation  for  the  integral  in 
(30)  is  accurate.  This  means  that  one  will  have  to  solve  a 
large  set  of  simultaneous  equations  in  order  to  produce  the 
desired  results.  The  solution  of  such  sets  of  equations  is 
a  machine  intensive  product  requiring  large  amounts  of  storage. 

At  this  time,  it  is  this  factor  which  restricts  the  major  possible 
numerical  solutions  of  the  Gelfand-Levitan  equation.  The  authors 
are  presently  investigating  the  possibility  of  lifting  or  at 
least  easing  this  restriction. 

Figures  12.  through  If  present  the  results  for  K(x,x)  versus  x  for 
four  potential  models:  a  single  delta  function  {12.  )y 
a  potential  with  a  rescaled  delta  function  reflection  coefficient  (1* 
a  2 -delta-function  potential  (13^  and  a  potential  step  (1 53 • 

The  results  exhibit  two  major  properties.  First  the  presence  of 
discontinuities  in  K(x,x)  does  not  impair  the  effectiveness  of  the 
algorithm,  as  the  discontinuities  one  expects  from  the  analytic 
solutions  for  V(x)  are  present  in  the  numerical  results.  Second, 
the  precision  of  the  results  gets  worse  as  x  increases.  This  is 
to  be  expected  since  any  quadrature  rule  used  to  approximate  the 
integral  in  (30)  will  be  less  accurate,  the  larger  the  range  of 
integration.  This  effect  can  be  reduced  by  keeping  the  grid  size 
small  and  using  those  approximations  which  are  best  suited  to  the 
form  of  the  integrand. 
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TABLE  I  -  IONOSPHERE  MODELS 


Model  I  -  Single  Delta  Function 

Potential 


Reflection 

Coefficient 


V (x)  =  Vq6(x) 

b(k)  =  -i((2k/Vo)+i)_1 


Fourier  transform  Figure  *7 
of  b (k) 


Virtual  height 

Gelfand-Levitan 

Solution 


Figure  J 
Figure  JZ 


Model  II  -  Rescaled  Delta  Function 

Potential  V(x)  =  vVq6(x)  -  VQa2sech2 (ox+a) 

o=  ( 1— v  2 ) 55  sinha  =  ((l-u)/2v)Js 

Reflection 

Coefficient  b(k)  =  -iv ( (2k/VQ) +i) _1 
Fourier  Transform 

of  b(k)  Figure  1  (rescale  by  v) 

Virtual  height  Figure  2 

Gelfand-Levitan 

Solution  Figure  14 


Model  III  -  Double  Delta  Function 

Potential  V(x)  =  Vq6  (x)  +  V^  (x-x^ 

Reflection 

Coefficient  b(k)  =  (- (1+ (2ik) /VQ) + (1- (2ik/Vi) ) e 

(-(l-(2ik/VQ) )  (l-(2ik)  A1)e“ 

Fourier  Transform 

of  b(k)  Figure  g 

Virtual  Height  Figure  £. 

Gelfand-Levitan 

Solution  Figure  13 


-2ikxi , 

x)  x 

2ikxi  +  1}-1 
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Model  IV  -  Potential  Step 

Potential  V(x)  =  n(x)V 


Reflection 

Coefficient 


U-i(B"2-l)3i)/(l+i(B"2- 


b(k) 


Fourier  Transform 
of  b (k) 

Virtual  Height 

Gelfand-Levitan 

Solution 


(l-(l-e''2)}s)/(l+(l-e'2) 


V'  6  +  k/ko 


Figure  9 
Figure  3 

Figure  15* 


% 


l)Js)6<  1 
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